Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db3) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Post-Transplant_1 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")


# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.1.4 Post-Transplant_2 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
          scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")
6.3.1.1.1 Number of samples used
[1] 27
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.2 Richness

betadisper(beta_div_richness_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.978
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.98
Hot_dry   0.97302        
adonis2(formula=beta_div_richness_wild$S ~ Population, data=wild[labels(beta_div_richness_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.542719 0.2095041 6.625717 0.001
Residual 25 5.820951 0.7904959 NA NA
Total 26 7.363669 1.0000000 NA NA

6.3.1.3 Neutral

betadisper(beta_div_neutral_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999   0.95
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.953
Hot_dry   0.94763        
adonis2(formula=beta_div_neutral_wild$S ~ Population, data=wild[labels(beta_div_neutral_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.918266 0.2608511 8.822682 0.001
Residual 25 5.435610 0.7391489 NA NA
Total 26 7.353876 1.0000000 NA NA

6.3.1.4 Phylogenetic

betadisper(beta_div_phylo_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.136
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.131
Hot_dry   0.12705        
adonis2(formula=beta_div_phylo_wild$S ~ Population, data=wild[labels(beta_div_phylo_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.3218613 0.2162815 6.899207 0.001
Residual 25 1.1662981 0.7837185 NA NA
Total 26 1.4881594 1.0000000 NA NA

6.3.1.5 Functional

betadisper(beta_div_func_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.019387 0.019387 1.653    999  0.199
Residuals 25 0.293200 0.011728                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.205
Hot_dry   0.21033        
adonis2(formula=beta_div_func_wild$S ~ Population, data=wild[labels(beta_div_func_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.0831048 0.1680538 5.05002 0.067
Residual 25 0.4114083 0.8319462 NA NA
Total 26 0.4945131 1.0000000 NA NA
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")
6.3.2.0.1 Number of samples used
[1] 26
beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.1 Richness

betadisper(beta_div_richness_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.093187 0.093187 11.812    999  0.003 **
Residuals 24 0.189340 0.007889                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.004
Hot_dry  0.0021532        
adonis2(formula=beta_div_richness_accli$S ~ Population, data=accli[labels(beta_div_richness_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.639630 0.1929088 5.736415 0.001
Residual 24 6.859879 0.8070912 NA NA
Total 25 8.499509 1.0000000 NA NA

6.3.2.2 Neutral

betadisper(beta_div_neutral_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05603 0.056026 4.1918    999   0.05 *
Residuals 24 0.32077 0.013365                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.048
Hot_dry  0.051717        
adonis2(formula=beta_div_neutral_accli$S ~ Population, data=accli[labels(beta_div_neutral_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 1.961192 0.2493171 7.970889 0.001
Residual 24 5.905063 0.7506829 NA NA
Total 25 7.866255 1.0000000 NA NA

6.3.2.3 Phylogenetic

betadisper(beta_div_phylo_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03637 0.036365 2.3087    999  0.149
Residuals 24 0.37804 0.015752                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.129
Hot_dry   0.14172        
adonis2(formula=beta_div_phylo_accli$S ~ Population, data=accli[labels(beta_div_phylo_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.2113395 0.1379515 3.84066 0.015
Residual 24 1.3206449 0.8620485 NA NA
Total 25 1.5319844 1.0000000 NA NA

6.3.2.4 Functional

betadisper(beta_div_func_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00409 0.004087 0.1789    999  0.681
Residuals 24 0.54821 0.022842                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.68
Hot_dry   0.67607        
adonis2(formula=beta_div_func_accli$S ~ Population, data=accli[labels(beta_div_func_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 1 0.0040177 0.0059732 0.1442177 0.63
Residual 24 0.6686040 0.9940268 NA NA
Total 25 0.6726217 1.0000000 NA NA
beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")
6.3.3.0.1 Number of samples used
[1] 53
beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
accli1_arrange<-accli1[labels(beta_div_neutral_accli1$S),]
accli1_arrange$Population_time <- interaction(accli1_arrange$Population, accli1_arrange$time_point)
6.3.3.0.2 Richness
betadisper(beta_div_richness_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.02724 0.0272402 3.7857    999  0.057 .
Residuals 51 0.36697 0.0071955                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                         0.065
1_Acclimation 0.057212              
adonis2(formula=beta_div_richness_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_richness_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.762933 0.2288365 4.846784 0.001
Residual 49 12.680830 0.7711635 NA NA
Total 52 16.443763 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 1.6396299 5.736415 0.19290877 0.001 0.006 *
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.5273732 1.906862 0.05462715 0.006 0.036 .
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 1.5558412 5.190707 0.17782052 0.001 0.006 *
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 1.8259388 8.319131 0.24968031 0.001 0.006 *
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.3856333 1.736034 0.09788177 0.009 0.054
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 1.5427188 6.625717 0.20950408 0.001 0.006 *
6.3.3.0.3 Neutral
betadisper(beta_div_neutral_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0083 0.0083036 0.9803    999   0.35
Residuals 51 0.4320 0.0084705                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              0_Wild 1_Acclimation
0_Wild                       0.339
1_Acclimation 0.3268              
adonis2(formula=beta_div_neutral_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_neutral_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 4.735485 0.2945657 6.820252 0.001
Residual 49 11.340673 0.7054343 NA NA
Total 52 16.076158 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 1.9611917 7.970889 0.24931708 0.001 0.006 *
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.7466371 3.064914 0.08498327 0.001 0.006 *
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 2.1421077 8.176866 0.25412251 0.001 0.006 *
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 2.0371666 10.078295 0.28730857 0.001 0.006 *
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.6314393 3.060027 0.16054681 0.001 0.006 *
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 1.9182663 8.822682 0.26085105 0.001 0.006 *
6.3.3.0.4 Phylogenetic
betadisper(beta_div_phylo_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00026 0.0002572 0.0158    999  0.903
Residuals 51 0.82986 0.0162718                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.904
1_Acclimation 0.90044              
adonis2(formula=beta_div_phylo_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_phylo_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.7622566 0.2345983 5.006223 0.001
Residual 49 2.4869430 0.7654017 NA NA
Total 52 3.2491997 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 0.2113395 3.840660 0.1379515 0.009 0.054
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.2076830 4.216762 0.1133028 0.003 0.018 .
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 0.3635093 4.799786 0.1666605 0.002 0.012 .
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 0.2121479 7.924059 0.2406769 0.001 0.006 *
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.2092433 3.885515 0.1953942 0.001 0.006 *
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 0.3218613 6.899207 0.2162815 0.001 0.006 *
6.3.3.0.5 Functional
betadisper(beta_div_func_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00006 0.0000604 0.0032    999  0.954
Residuals 51 0.95116 0.0186502                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.954
1_Acclimation 0.95485              
adonis2(formula=beta_div_func_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_func_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.1024301 0.0866259 1.549079 0.338
Residual 49 1.0800123 0.9133741 NA NA
Total 52 1.1824424 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 0.0040176901 0.144217742 0.0059731793 0.646 1.000
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.0002775823 0.009969841 0.0003020251 0.698 1.000
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 0.0871206378 3.762920746 0.1355376396 0.067 0.402
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 0.0025235151 0.120315257 0.0047895600 0.595 1.000
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.0409732757 4.066330477 0.2026444487 0.127 0.762
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 0.0831047968 5.050019518 0.1680537849 0.048 0.288
beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Effect of FMT on microbiota community

6.3.4.1 Comparison between Acclimation vs Post-FMT1

post3 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post3))]
identical(sort(colnames(post3.counts)),sort(as.character(rownames(post3))))

post3_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")
6.3.4.1.1 Number of samples used
[1] 52
beta_div_richness_post3<-hillpair(data=post3.counts, q=0)
beta_div_neutral_post3<-hillpair(data=post3.counts, q=1)
beta_div_phylo_post3<-hillpair(data=post3.counts, q=1, tree=genome_tree)
beta_div_func_post3<-hillpair(data=post3.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post3_arrange<-post3[labels(beta_div_neutral_post3$S),]
post3_arrange$type_time <- interaction(post3_arrange$type, post3_arrange$time_point)
6.3.4.1.2 Richness
betadisper(beta_div_richness_post3$S, post3_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     2 0.087391 0.043695 9.5886    999  0.001 ***
Residuals 49 0.223293 0.004557                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Control Hot_control Treatment
Control                 0.00400000     0.524
Hot_control 0.00219320                 0.001
Treatment   0.49723207  0.00010148          
adonis2(formula=beta_div_richness_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 4.399457 0.2450497 2.986233 0.001
Residual 46 13.553865 0.7549503 NA NA
Total 51 17.953321 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_richness_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.224 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.6804988 2.100611 0.12283837 0.002 0.030 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8481258 2.673290 0.16033371 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.1180661 3.809275 0.20252111 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.7216200 2.172734 0.11956009 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9551308 2.926054 0.16322910 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.2263345 4.039487 0.20157637 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.4319792 5.384836 0.25180628 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8172413 3.194690 0.17558364 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.5796135 2.441615 0.13239702 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.010 0.150
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.106 1.000
6.3.4.1.3 Neutral
betadisper(beta_div_neutral_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00430 0.0043011 0.5537    999  0.473
Residuals 50 0.38839 0.0077677                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                      0.47
5_Post-FMT1         0.46029            
adonis2(formula=beta_div_neutral_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 5.32464 0.3149916 4.230493 0.001
Residual 46 11.57943 0.6850084 NA NA
Total 51 16.90407 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_neutral_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.460 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.8153894 3.0970603 0.17113610 0.001 0.015 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.1809241 4.4856470 0.24265567 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.4321524 5.7774260 0.27806264 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.9517634 3.3715700 0.17404733 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 1.3127773 4.6298256 0.23585668 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.6713369 6.2395460 0.28056085 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.5409781 6.8338056 0.29928456 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9133614 4.0964534 0.21451383 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.6954835 3.2951234 0.17077493 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.2508491 0.13047758 0.014 0.210
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.1436369 0.20570451 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.6372683 0.09840968 0.040 0.600
6.3.4.1.4 Phylogenetic
betadisper(beta_div_phylo_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02481 0.024807 1.6411    999  0.236
Residuals 50 0.75581 0.015116                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.226
5_Post-FMT1         0.20609            
adonis2(formula=beta_div_phylo_post3$S ~ time_point*type, data=post3[labels(beta_div_phylo_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.7221892 0.2434166 2.959929 0.002
Residual 46 2.2446957 0.7565834 NA NA
Total 51 2.9668849 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_phylo_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.7893300 0.04999136 0.540 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.0547314 0.16919285 0.008 0.120
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.11760287 2.9988032 0.16661126 0.020 0.300
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09624422 1.7247115 0.10968160 0.176 1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.18586788 4.1904227 0.21836011 0.005 0.075
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.004 0.060
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.26358465 4.3608960 0.21417997 0.003 0.045 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.25319427 3.2738422 0.17915456 0.045 0.675
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.39050120 5.9837393 0.27218933 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 0.14203376 5.4200212 0.25303529 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09666753 2.3682173 0.13635351 0.013 0.195
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.09252600 2.9824958 0.15711821 0.005 0.075
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.784 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.108 1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.694 1.000
6.3.4.1.5 Functional
betadisper(beta_div_func_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02408 0.024081 1.2619    999  0.251
Residuals 50 0.95415 0.019083                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.245
5_Post-FMT1         0.26666            
adonis2(formula=beta_div_func_post3$S ~ time_point*type, data=post3[labels(beta_div_func_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.1702325 0.1167018 1.215508 0.245
Residual 46 1.2884648 0.8832982 NA NA
Total 51 1.4586973 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_func_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.022076664 0.65002192 0.041534889 0.430 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.016749605 0.52091146 0.033561912 0.479 1
Control.1_Acclimation vs Control.5_Post-FMT1 1 -0.008325555 -0.22800110 -0.015434681 0.900 1
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.126527918 3.45550519 0.197960767 0.076 1
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.051008429 1.23951838 0.076327287 0.283 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.001022590 0.05430389 0.003382513 0.662 1
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.002157067 0.09411569 0.005847832 0.618 1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.056602363 2.56037069 0.145803909 0.179 1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.009569124 0.35095521 0.021463896 0.495 1
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 -0.001745663 -0.08225018 -0.005167199 0.703 1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.057758674 2.84545622 0.159449901 0.148 1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.005575266 0.21803560 0.013444020 0.559 1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.119540855 4.84764704 0.244242909 0.079 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.052587837 1.77308932 0.099762584 0.215 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.012980354 0.44307662 0.028690955 0.466 1
beta_richness_nmds_post3 <- beta_div_richness_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- beta_div_neutral_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- beta_div_phylo_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- beta_div_func_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

6.3.4.2 Comparison between Acclimation vs Post-FMT2

post4 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post4.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post4))]
identical(sort(colnames(post4.counts)),sort(as.character(rownames(post4))))

post4_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")
6.3.4.2.1 Number of samples used
[1] 53
beta_div_richness_post4<-hillpair(data=post4.counts, q=0)
beta_div_neutral_post4<-hillpair(data=post4.counts, q=1)
beta_div_phylo_post4<-hillpair(data=post4.counts, q=1, tree=genome_tree)
beta_div_func_post4<-hillpair(data=post4.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post4_arrange<-post4[labels(beta_div_neutral_post4$S),]
post4_arrange$type_time <- interaction(post4_arrange$type, post4_arrange$time_point)
6.3.4.2.2 Richness
betadisper(beta_div_richness_post4$S, post4_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     2 0.05757 0.0287844 3.754    999  0.031 *
Residuals 50 0.38338 0.0076677                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0840000     0.497
Hot_control 0.0769856                 0.003
Treatment   0.4949731   0.0047188          
adonis2(formula=beta_div_richness_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.279021 0.1937289 3.924535 0.001
Residual 49 13.646799 0.8062711 NA NA
Total 52 16.925820 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_richness_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.251 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.8072604 2.940901 0.16392161 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4817387 1.660775 0.09968176 0.027 0.405
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.1179704 3.885459 0.20573812 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.9130048 3.195028 0.16645080 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.5959230 1.984036 0.11032208 0.004 0.060
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.2747787 4.275366 0.21086503 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6397330 2.913695 0.15405213 0.002 0.030 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.4575447 6.224524 0.28007456 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3276169 1.412318 0.08111028 0.034 0.510
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.4.2.3 Neutral
betadisper(beta_div_neutral_post4$S, post4_arrange$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00848 0.0084800 1.2368    999  0.269
Residuals 51 0.34967 0.0068563                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                      0.26
6_Post-FMT2          0.2713            
adonis2(formula=beta_div_neutral_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 3.823089 0.2357666 5.038847 0.001
Residual 49 12.392477 0.7642334 NA NA
Total 52 16.215567 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_neutral_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.479 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 1.1524353 4.9536068 0.24825621 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4748999 1.9609749 0.11561688 0.006 0.090
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.3503168 5.4081420 0.26499923 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 1.3540292 5.3398081 0.25022756 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.6311089 2.4041625 0.13063146 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.6125755 5.9825981 0.27215155 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6202327 3.1519868 0.16457754 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.5701179 7.6327037 0.32297209 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3634438 1.7083388 0.09647087 0.030 0.450
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.6483346 0.22511910 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.2065321 0.12119453 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.7710313 0.26507845 0.001 0.015 .
6.3.4.2.4 Phylogenetic
betadisper(beta_div_phylo_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03769 0.037692 3.8614    999  0.045 *
Residuals 51 0.49783 0.009761                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                      0.04
6_Post-FMT2        0.054869            
adonis2(formula=beta_div_phylo_post4$S ~ time_point*Population, data=post4[labels(beta_div_phylo_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.6466392 0.2442162 5.277784 0.001
Residual 49 2.0011759 0.7557838 NA NA
Total 52 2.6478151 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_phylo_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.789330 0.04999136 0.548 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.054731 0.16919285 0.015 0.225
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.15591805 4.379209 0.22597458 0.005 0.075
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.07099742 1.879364 0.11134092 0.124 1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.18682367 4.878754 0.24542555 0.003 0.045 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.052184 0.20208192 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.36496892 6.396667 0.28560797 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.22628210 3.829222 0.19311005 0.015 0.225
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.34830814 5.846334 0.26761166 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.10002871 4.383624 0.21505615 0.002 0.030 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 0.12577510 5.060129 0.24027055 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.06334378 2.499774 0.13512455 0.021 0.315
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.382025 0.12958449 0.028 0.420
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.722460 0.14541146 0.007 0.105
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.043656 0.20174244 0.001 0.015 .
6.3.4.2.5 Functional
betadisper(beta_div_func_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00024 0.0002438 0.0119    999  0.925
Residuals 51 1.04336 0.0204581                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.926
6_Post-FMT2         0.91349            
adonis2(formula=beta_div_func_post4$S ~ time_point*Population, data=post4[labels(beta_div_func_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 3 0.0258537 0.0188483 0.313769 0.578
Residual 49 1.3458235 0.9811517 NA NA
Total 52 1.3716772 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_func_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.0220766636 0.650021923 0.0415348890 0.407 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.0167496054 0.520911459 0.0335619117 0.478 1
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.0376809680 1.119650309 0.0694587220 0.307 1
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.0292200956 0.920511083 0.0578191917 0.390 1
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0133433458 0.270473784 0.0177122064 0.510 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.0010225901 0.054303886 0.0033825127 0.630 1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 -0.0005177706 -0.025585400 -0.0016016487 0.727 1
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.0013301207 0.072110871 0.0044867082 0.600 1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0060959077 0.174487757 0.0107878382 0.613 1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.0010345754 0.055797964 0.0034752533 0.628 1
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 -0.0001056284 -0.006306177 -0.0003942915 0.689 1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0017235602 0.051851181 0.0032302306 0.755 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -0.0080428882 -0.442986255 -0.0284750185 0.851 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -0.0011796256 -0.034047378 -0.0021324990 0.904 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.0036300838 0.110487573 0.0068581148 0.712 1
beta_richness_nmds_post4 <- beta_div_richness_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- beta_div_neutral_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- beta_div_phylo_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- beta_div_func_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

6.3.4.3 Comparison between Acclimation vs Post-FMT1 and Post-FMT2

post6 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post6.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post6))]
identical(sort(colnames(post6.counts)),sort(as.character(rownames(post6))))

post6_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")
6.3.4.3.1 Number of samples used
[1] 79
beta_div_richness_post6<-hillpair(data=post6.counts, q=0)
beta_div_neutral_post6<-hillpair(data=post6.counts, q=1)
beta_div_phylo_post6<-hillpair(data=post6.counts, q=1, tree=genome_tree)
beta_div_func_post6<-hillpair(data=post6.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post6_arrange<-post6[labels(beta_div_neutral_post6$S),]
post6_arrange$type_time <- interaction(post6_arrange$type, post6_arrange$time_point)
6.3.4.3.2 Richness
betadisper(beta_div_richness_post6$S, post6_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06454 0.032268 4.9753    999  0.014 *
Residuals 76 0.49290 0.006485                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0120000     0.681
Hot_control 0.0119513                 0.003
Treatment   0.6817089   0.0035106          
adonis2(formula=beta_div_richness_post6$S ~ time_point*type, data=post6[labels(beta_div_richness_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 8 6.395467 0.2451322 2.841433 0.001
Residual 70 19.694403 0.7548678 NA NA
Total 78 26.089870 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_richness_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.244 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.036 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.6804988 2.100611 0.12283837 0.001 0.036 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8481258 2.673290 0.16033371 0.001 0.036 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.1180661 3.809275 0.20252111 0.001 0.036 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.8072604 2.940901 0.16392161 0.001 0.036 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4817387 1.660775 0.09968176 0.029 1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.1179704 3.885459 0.20573812 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.036 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.7216200 2.172734 0.11956009 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9551308 2.926054 0.16322910 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.2263345 4.039487 0.20157637 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.9130048 3.195028 0.16645080 0.002 0.072
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.5959230 1.984036 0.11032208 0.002 0.072
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.2747787 4.275366 0.21086503 0.002 0.072
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.4319792 5.384836 0.25180628 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8172413 3.194690 0.17558364 0.002 0.072
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.5796135 2.441615 0.13239702 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6397330 2.913695 0.15405213 0.001 0.036 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.4575447 6.224524 0.28007456 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3276169 1.412318 0.08111028 0.036 1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.011 0.396
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.036 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.7628135 2.683925 0.14364890 0.001 0.036 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3432605 1.148733 0.06698647 0.266 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.1269580 3.799256 0.19188884 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.110 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3571397 1.297184 0.07959561 0.122 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.7769467 2.670898 0.15114670 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.6502360 2.253407 0.13060650 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4132091 1.616138 0.09174188 0.008 0.288
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0163992 3.760571 0.19030682 0.002 0.072
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.2732563 1.019281 0.05988979 0.466 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.002 0.072
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.001 0.036 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.036 .
6.3.4.3.3 Neutral
betadisper(beta_div_neutral_post6$S, post6_arrange$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00896 0.0044799 0.6329    999  0.508
Residuals 76 0.53799 0.0070788                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                   0.45700       0.251
5_Post-FMT1         0.45960                   0.756
6_Post-FMT2         0.27047     0.75275            
adonis2(formula=beta_div_neutral_post6$S ~ time_point*type, data=post6[labels(beta_div_neutral_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 8 7.736408 0.3121974 3.971673 0.001
Residual 70 17.044094 0.6878026 NA NA
Total 78 24.780501 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_neutral_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.453 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.036 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.8153894 3.0970603 0.17113610 0.001 0.036 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.1809241 4.4856470 0.24265567 0.001 0.036 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.4321524 5.7774260 0.27806264 0.001 0.036 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 1.1524353 4.9536068 0.24825621 0.001 0.036 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4748999 1.9609749 0.11561688 0.004 0.144
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.3503168 5.4081420 0.26499923 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.036 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.9517634 3.3715700 0.17404733 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 1.3127773 4.6298256 0.23585668 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.6713369 6.2395460 0.28056085 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 1.3540292 5.3398081 0.25022756 0.001 0.036 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.6311089 2.4041625 0.13063146 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.6125755 5.9825981 0.27215155 0.001 0.036 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.5409781 6.8338056 0.29928456 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9133614 4.0964534 0.21451383 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.6954835 3.2951234 0.17077493 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6202327 3.1519868 0.16457754 0.001 0.036 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.5701179 7.6327037 0.32297209 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3634438 1.7083388 0.09647087 0.041 1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.2508491 0.13047758 0.011 0.396
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.1436369 0.20570451 0.001 0.036 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.8908158 3.7146920 0.18842252 0.002 0.072
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3860927 1.5521758 0.08843210 0.085 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.3122237 5.1302726 0.24279254 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.6372683 0.09840968 0.045 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3157079 1.3252026 0.08117526 0.152 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0579520 4.2700097 0.22158835 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.7454015 2.9200493 0.16294873 0.002 0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4377161 1.9421261 0.10824392 0.009 0.324
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.3766597 5.8752789 0.26858075 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.3176516 1.3161367 0.07600637 0.169 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.6483346 0.22511910 0.001 0.036 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.2065321 0.12119453 0.003 0.108
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.7710313 0.26507845 0.001 0.036 .
6.3.4.3.4 Phylogenetic
betadisper(beta_div_phylo_post6$S, post6$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.04250 0.021252 1.8628    999  0.161
Residuals 76 0.86703 0.011408                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                  0.223000       0.051
5_Post-FMT1        0.205952                   0.746
6_Post-FMT2        0.054864    0.717296            
adonis2(formula=beta_div_phylo_post6$S ~ time_point*type, data=post6[labels(beta_div_phylo_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 8 1.055866 0.2692287 3.223651 0.001
Residual 70 2.865952 0.7307713 NA NA
Total 78 3.921818 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_phylo_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.7893300 0.04999136 0.545 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.0547314 0.16919285 0.010 0.360
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.11760287 2.9988032 0.16661126 0.024 0.864
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09624422 1.7247115 0.10968160 0.163 1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.18586788 4.1904227 0.21836011 0.003 0.108
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.15591805 4.3792085 0.22597458 0.005 0.180
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.07099742 1.8793639 0.11134092 0.120 1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.18682367 4.8787543 0.24542555 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.003 0.108
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.26358465 4.3608960 0.21417997 0.002 0.072
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.25319427 3.2738422 0.17915456 0.044 1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.39050120 5.9837393 0.27218933 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.36496892 6.3966666 0.28560797 0.001 0.036 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.22628210 3.8292220 0.19311005 0.015 0.540
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.34830814 5.8463335 0.26761166 0.002 0.072
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 0.14203376 5.4200212 0.25303529 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09666753 2.3682173 0.13635351 0.014 0.504
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.09252600 2.9824958 0.15711821 0.012 0.432
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.10002871 4.3836237 0.21505615 0.002 0.072
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 0.12577510 5.0601287 0.24027055 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.06334378 2.4997737 0.13512455 0.014 0.504
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.795 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.102 1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07917244 3.0180046 0.15869197 0.006 0.216
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.04335491 1.5335604 0.08746429 0.162 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.10783045 3.7500438 0.18987521 0.002 0.072
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.704 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.06393539 1.5651817 0.09448624 0.132 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.05265949 1.2240203 0.07544494 0.295 1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.09753501 2.2402429 0.12994265 0.019 0.684
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07228545 2.3279593 0.12701683 0.035 1.000
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.11759094 3.5538444 0.18174658 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.06667255 1.9859527 0.11041687 0.089 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.028 1.000
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.004 0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.002 0.072
6.3.4.3.5 Functional
betadisper(beta_div_func_post6$S, post6$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.03015 0.015074 0.7893    999  0.487
Residuals 76 1.45147 0.019098                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                   0.29300       0.932
5_Post-FMT1         0.26551                   0.298
6_Post-FMT2         0.93087     0.28011            
adonis2(formula=beta_div_func_post6$S ~ time_point*type, data=post6[labels(beta_div_func_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()
Df SumOfSqs R2 F Pr(>F)
Model 8 0.1632318 0.0763851 0.7236452 0.472
Residual 70 1.9737271 0.9236149 NA NA
Total 78 2.1369589 1.0000000 NA NA
pairwise <- pairwise.adonis(beta_div_func_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.0220766636 0.650021923 0.0415348890 0.456 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.0167496054 0.520911459 0.0335619117 0.483 1
Control.1_Acclimation vs Control.5_Post-FMT1 1 -0.0083255553 -0.228001097 -0.0154346814 0.882 1
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.1265279179 3.455505193 0.1979607668 0.077 1
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.0510084290 1.239518381 0.0763272870 0.297 1
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.0376809680 1.119650309 0.0694587220 0.304 1
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.0292200956 0.920511083 0.0578191917 0.374 1
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0133433458 0.270473784 0.0177122064 0.505 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.0010225901 0.054303886 0.0033825127 0.645 1
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.0021570675 0.094115687 0.0058478321 0.635 1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.0566023632 2.560370692 0.1458039091 0.161 1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.0095691236 0.350955208 0.0214638964 0.509 1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 -0.0005177706 -0.025585400 -0.0016016487 0.715 1
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.0013301207 0.072110871 0.0044867082 0.592 1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0060959077 0.174487757 0.0107878382 0.596 1
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 -0.0017456629 -0.082250179 -0.0051671989 0.719 1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.0577586745 2.845456220 0.1594499006 0.146 1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.0055752661 0.218035597 0.0134440201 0.548 1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.0010345754 0.055797964 0.0034752533 0.659 1
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 -0.0001056284 -0.006306177 -0.0003942915 0.673 1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.0017235602 0.051851181 0.0032302306 0.741 1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.1195408549 4.847647043 0.2442429086 0.072 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0525878365 1.773089316 0.0997625840 0.226 1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0265995825 1.175418056 0.0684360667 0.297 1
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0145818992 0.699759916 0.0419023938 0.414 1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0080695208 -0.216173226 -0.0136958691 0.931 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0129803540 0.443076619 0.0286909552 0.461 1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0267162134 1.225605811 0.0755352882 0.272 1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0384388433 1.932815819 0.1141461550 0.241 1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.0553988290 1.478193905 0.0897060633 0.230 1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 -0.0040061386 -0.148504693 -0.0093684974 0.744 1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0024023972 0.095389804 0.0059265296 0.629 1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0004960759 -0.011903277 -0.0007445087 0.849 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -0.0080428882 -0.442986255 -0.0284750185 0.860 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -0.0011796256 -0.034047378 -0.0021324990 0.917 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.0036300838 0.110487573 0.0068581148 0.686 1
beta_richness_nmds_post6 <- beta_div_richness_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post6 <- beta_div_neutral_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post6 <- beta_div_phylo_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post6 <- beta_div_func_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Are there differences between the control and the treatment group?

6.3.5.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")
6.3.5.1.1 Number of samples used
[1] 26
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]
6.3.5.1.2 Richness
betadisper(beta_div_richness_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.017675 0.0088373 2.3825    999  0.097 .
Residuals 23 0.085312 0.0037092                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0040000     0.647
Hot_control 0.0068795                 0.231
Treatment   0.6248469   0.2084296          
adonis2(formula=beta_div_richness_post1$S ~ type, data=post1[labels(beta_div_richness_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.195567 0.1448246 1.947534 0.001
Residual 23 7.059710 0.8551754 NA NA
Total 25 8.255277 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.010      0.030   .
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.108      0.324    
6.3.5.1.3 Neutral
betadisper(beta_div_neutral_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999   0.58
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.20200     0.957
Hot_control 0.21166                 0.445
Treatment   0.95468     0.43604          
adonis2(formula=beta_div_neutral_post1$S ~ type, data=post1[labels(beta_div_neutral_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.395968 0.1900228 2.697931 0.001
Residual 23 5.950350 0.8099772 NA NA
Total 25 7.346318 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.012      0.036   .
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.001      0.003   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.061      0.183    
6.3.5.1.4 Phylogenetic
betadisper(beta_div_phylo_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.911
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92200     0.692
Hot_control 0.91505                 0.770
Treatment   0.63312     0.73046          
adonis2(formula=beta_div_phylo_post1$S ~ type, data=post1[labels(beta_div_phylo_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.0745104 0.0705947 0.8735033 0.537
Residual 23 0.9809570 0.9294053 NA NA
Total 25 1.0554673 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.771       1.00    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.110       0.33    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.694       1.00    
6.3.5.1.5 Functional
betadisper(beta_div_func_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00400 0.0019999 0.1431    999  0.859
Residuals 23 0.32135 0.0139717                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.59500     0.744
Hot_control 0.60188                 0.813
Treatment   0.74597     0.84473          
adonis2(formula=beta_div_func_post1$S ~ type, data=post1[labels(beta_div_func_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.1230554 0.1608583 2.204479 0.174
Residual 23 0.6419374 0.8391417 NA NA
Total 25 0.7649929 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post1$S,post1_arrange$type, perm=999)
pairwise
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11954085 4.8476470 0.24424291   0.068      0.204    
2   Control vs Hot_control  1 0.05258784 1.7730893 0.09976258   0.227      0.681    
3 Treatment vs Hot_control  1 0.01298035 0.4430766 0.02869096   0.453      1.000    
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")
6.3.5.2.1 Number of samples used
[1] 27
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]
6.3.5.2.2 Richness
betadisper(beta_div_richness_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.845
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.72000     0.818
Hot_control 0.67789                 0.636
Treatment   0.79246     0.59820          
adonis2(formula=beta_div_richness_post2$S ~ type, data=post2[labels(beta_div_richness_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.504341 0.1967776 2.939822 0.001
Residual 24 6.140538 0.8032224 NA NA
Total 26 7.644879 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.6463814 2.560441 0.1379515 0.001 0.003 *
Treatment vs Hot_control 1 0.4796256 1.916520 0.1069694 0.001 0.003 *
Control vs Hot_control 1 1.1305044 4.268317 0.2105906 0.001 0.003 *
6.3.5.2.3 Neutral
betadisper(beta_div_neutral_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.434
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.43400     0.651
Hot_control 0.44675                 0.236
Treatment   0.65989     0.25095          
adonis2(formula=beta_div_neutral_post2$S ~ type, data=post2[labels(beta_div_neutral_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 1.923807 0.2603795 4.224537 0.001
Residual 24 5.464666 0.7396205 NA NA
Total 26 7.388473 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 1.0227481 4.648335 0.2251191 0.001 0.003 *
Treatment vs Hot_control 1 0.5010202 2.206532 0.1211945 0.001 0.003 *
Control vs Hot_control 1 1.3619424 5.771031 0.2650785 0.001 0.003 *
6.3.5.2.4 Phylogenetic
betadisper(beta_div_phylo_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.956
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93900     0.834
Hot_control 0.93765                 0.771
Treatment   0.83933     0.76015          
adonis2(formula=beta_div_phylo_post2$S ~ type, data=post2[labels(beta_div_phylo_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 0.1594363 0.2042241 3.079623 0.001
Residual 24 0.6212564 0.7957759 NA NA
Total 26 0.7806927 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.05927454 2.382025 0.1295845 0.033 0.099
Treatment vs Hot_control 1 0.06906280 2.722460 0.1454115 0.003 0.009 *
Control vs Hot_control 1 0.11081709 4.043656 0.2017424 0.003 0.009 *
6.3.5.2.5 Functional
betadisper(beta_div_func_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01259 0.0062962 0.3249    999  0.777
Residuals 24 0.46507 0.0193778                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.50600     0.646
Hot_control 0.45381                 0.776
Treatment   0.57452     0.74365          
adonis2(formula=beta_div_func_post2$S ~ type, data=post2[labels(beta_div_func_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 2 -0.0037283 -0.0054704 -0.065288 0.923
Residual 24 0.6852623 1.0054704 NA NA
Total 26 0.6815340 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 -0.008042888 -0.44298625 -0.028475019 0.847 1
Treatment vs Hot_control 1 -0.001179626 -0.03404738 -0.002132499 0.894 1
Control vs Hot_control 1 0.003630084 0.11048757 0.006858115 0.684 1
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")
6.3.5.3.1 Number of samples used
[1] 53
beta_div_richness_post5<-hillpair(data=post5.counts, q=0)
beta_div_neutral_post5<-hillpair(data=post5.counts, q=1)
beta_div_phylo_post5<-hillpair(data=post5.counts, q=1, tree=genome_tree)
beta_div_func_post5<-hillpair(data=post5.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post5_arrange<-post5[labels(beta_div_neutral_post5$S),]
post5_arrange$type_time <- interaction(post5_arrange$type, post5_arrange$time_point)
6.3.5.3.2 Richness
betadisper(beta_div_richness_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999  0.186
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.041000     0.715
Hot_control 0.039117                 0.228
Treatment   0.716358    0.218648          
adonis2(formula=beta_div_richness_post5$S ~ type*time_point, data=post5[labels(beta_div_richness_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 3.261916 0.1981463 2.322836 0.005
Residual 47 13.200248 0.8018537 NA NA
Total 52 16.462165 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_richness_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.014 0.210
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.7628135 2.683925 0.14364890 0.001 0.015 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3432605 1.148733 0.06698647 0.270 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.1269580 3.799256 0.19188884 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.111 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3571397 1.297184 0.07959561 0.137 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.7769467 2.670898 0.15114670 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.6502360 2.253407 0.13060650 0.003 0.045 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4132091 1.616138 0.09174188 0.003 0.045 .
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0163992 3.760571 0.19030682 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.2732563 1.019281 0.05988979 0.423 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.5.3.3 Neutral
betadisper(beta_div_neutral_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999   0.23
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.10800     0.884
Hot_control 0.10701                 0.178
Treatment   0.87156     0.17449          
adonis2(formula=beta_div_neutral_post5$S ~ type*time_point, data=post5[labels(beta_div_neutral_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 3.947979 0.2569798 3.251069 0.002
Residual 47 11.415016 0.7430202 NA NA
Total 52 15.362995 1.0000000 NA NA
pairwise<-pairwise.adonis(beta_div_neutral_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.250849 0.13047758 0.020 0.300
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.143637 0.20570451 0.002 0.030 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.8908158 3.714692 0.18842252 0.001 0.015 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3860927 1.552176 0.08843210 0.097 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.3122237 5.130273 0.24279254 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.637268 0.09840968 0.043 0.645
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3157079 1.325203 0.08117526 0.156 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0579520 4.270010 0.22158835 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.7454015 2.920049 0.16294873 0.002 0.030 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4377161 1.942126 0.10824392 0.006 0.090
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.3766597 5.875279 0.26858075 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.3176516 1.316137 0.07600637 0.181 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.648335 0.22511910 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.206532 0.12119453 0.002 0.030 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.771031 0.26507845 0.001 0.015 .
6.3.5.3.4 Phylogenetic
betadisper(beta_div_phylo_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999  0.985
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.89500     0.855
Hot_control 0.88926                 0.930
Treatment   0.82391     0.91902          
adonis2(formula=beta_div_phylo_post5$S ~ type*time_point, data=post5[labels(beta_div_phylo_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.3518222 0.180049 2.0641 0.046
Residual 47 1.6022134 0.819951 NA NA
Total 52 1.9540356 1.000000 NA NA
pairwise<-pairwise.adonis(beta_div_phylo_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.769 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.108 1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07917244 3.0180046 0.15869197 0.005 0.075
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.04335491 1.5335604 0.08746429 0.193 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.10783045 3.7500438 0.18987521 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.698 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.06393539 1.5651817 0.09448624 0.173 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.05265949 1.2240203 0.07544494 0.299 1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.09753501 2.2402429 0.12994265 0.010 0.150
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07228545 2.3279593 0.12701683 0.037 0.555
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.11759094 3.5538444 0.18174658 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.06667255 1.9859527 0.11041687 0.101 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.024 0.360
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.003 0.045 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.003 0.045 .
6.3.5.3.5 Functional
betadisper(beta_div_func_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00785 0.0039232 0.2322    999  0.811
Residuals 50 0.84483 0.0168966                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.54500     0.623
Hot_control 0.52384                 0.835
Treatment   0.58787     0.85068          
adonis2(formula=beta_div_func_post5$S ~ type*time_point, data=post5[labels(beta_div_func_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()
Df SumOfSqs R2 F Pr(>F)
Model 5 0.1076682 0.075037 0.7625685 0.465
Residual 47 1.3271997 0.924963 NA NA
Total 52 1.4348679 1.000000 NA NA
pairwise<-pairwise.adonis(beta_div_func_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.1195408549 4.84764704 0.2442429086 0.067 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0525878365 1.77308932 0.0997625840 0.264 1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0265995825 1.17541806 0.0684360667 0.312 1
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0145818992 0.69975992 0.0419023938 0.388 1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0080695208 -0.21617323 -0.0136958691 0.923 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.0129803540 0.44307662 0.0286909552 0.456 1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.0267162134 1.22560581 0.0755352882 0.337 1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0384388433 1.93281582 0.1141461550 0.218 1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.0553988290 1.47819391 0.0897060633 0.250 1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 -0.0040061386 -0.14850469 -0.0093684974 0.729 1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.0024023972 0.09538980 0.0059265296 0.612 1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -0.0004960759 -0.01190328 -0.0007445087 0.842 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -0.0080428882 -0.44298625 -0.0284750185 0.838 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -0.0011796256 -0.03404738 -0.0021324990 0.912 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.0036300838 0.11048757 0.0068581148 0.681 1
beta_richness_nmds_post5 <- beta_div_richness_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- beta_div_neutral_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- beta_div_phylo_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- beta_div_func_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")